Surface Integral Method for the Second Harmonic Generation in Metal Nanoparticles 
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Second harmonic (SH) radiation in metal nanoparticles is generated by both nonlocal-bulk and 
local-surface SH sources, induced by the electromagnetic field at the fundamental frequency. We 
propose a surface integral equation (SIE) method for evaluating the SH radiation generated by metal 
nanoparticles with arbitrary shapes, considering both the local-surface and nonlocal-bulk SH sources. 
We demonstrate that the contribution of the nonlocal-bulk SH source to the SH electromagnetic field 
can be rigorously taken into account through an equivalent surface magnetic current. We numerically 
solve the SIE problem by using the Galerkin method and the Rao- Wilton- Glisson basis functions 
in the framework of the distribution theory. The accuracy of the proposed method is verified by 
comparison with the SH-Mie formulation in the case of SH generation from a spherical nanoparticle. 
As an example of a complex-shaped particle, we investigate the SH radiation generated from a 
triangular nano-prism. The method paves the way for a better understanding of the SH generation 
process in arbitrarily shaped nanoparticles and can also have a high impact in the design of novel 
nanoplasmonic devices with enhanced SH emission. 
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I. INTRODUCTION 

The study of nonlinear optical properties of noble 
metal nanoparticles (NPs) has become a very active re- 
search field, stimulated by both the maturity of nanofab- 
rication techniques and the appeal of the potential ap- 
plications. Noble metal nanoparticles support localized 
surface plasmon (LSP) resonances, consisting in the col- 
lective oscillations of the conduction electrons [l[. LSPs 
can significantly enhance both the linear and nonlinear 
scattering processes @ , enabling strong nonlinear optical 
processes at relatively low excitation powers. 

The second harmonic (SH) generation is a nonlinear 
optical process in which two photons of the same fre- 
quency, i.e. the fundamental frequency, interact with 
the material and generate one photon of twice the fre- 
quency. SH generation from noble metal NPs takes 
place in the bulk of the particle and on its surface. In 
particular, the electromagnetic field at the fundamen- 
tal frequency induces nonlocal-bulk and local-surface SH 
sources, whereas the local-bulk SH source is absent due 
to the centrosymmetry of the material [3J. The local- 
surface SH source is due to the symmetry breaking at the 
interface between the NP and the embedding medium Q ■ 
The relative contribution to the SH radiation of nonlocal- 
bulk and local-surface sources depend on the NP shape 
and size, and on the optical properties of the metal at 
both the fundamental and the SH frequencies [5|. 

Several analytical models have been developed to de- 
scribe the SHG from noble metals NPs of simple shapes 
such as spheres and cylinders. In particular, Dadap et. 
al [|| studied the SH radiation generated from the local- 
surface SH source in a sphere with radius much smaller 
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than the wavelength of the exciting light (Rayleigh limit). 
This model was later extended in ja0|) considering also 
the nonlocal-bulk source. A full-wave theory of the SH 
generation by a sphere of centrosymmetric material fea- 
turing SH local-surface sources has been proposed in [1] . 
This approach has been recently extended to model ag- 
gregates of spheres @ but, also in this case, only the 
local-surface SH sources have been taken into account. 
A full-wave analytical theory has been also proposed for 
the SH scattering of light by a cylinder of infinite length 
10]. 

The last few decades have witnessed a continuous im- 
provement of nanofabrication techniques and a large vari- 
ety of complex shapes can be controllably achieved nowa- 
days. This fact allows for an accurate tuning of the spec- 
tral position of the LSP resonance of a NP and stimulates 
the demand for numerical tools to model the linear [ll[ 
and non-linear behavior of NPs with arbitrary shapes. In 
fact, De Beer et al. [H[ presented a theoretical frame- 
work for solving the SH generation problem for particles 
with arbitrary shapes in the low contrast limit. Makitalo 
et. al. [HI introduced a surface integral equation (SIE) 
formulation for studying the nonlinear scattering from 
arbitrarily shaped particles taking into account only the 
local-surface SH sources. However, the validity of their 
numerical formulation will be questioned in the present 
paper. Benedetti et. al. [HI developed a volume integral 
equation (VIE) formulation based on the dyadic Green 
function for the study of non-linear scattering properties 
from arbitrarily shaped particles due to both nonlocal- 
bulk and local-surface contributions. However, unlike 
SIE formulations, the VIE formulations require the dis- 
cretization of the entire NP volume, resulting in an higher 
memory occupancy and longer computational time. 

In this paper, we present a SIE formulation of the 
SH generation problem for three dimensional, arbitrar- 
ily shaped nanoparticles made of centrosymmetric and 
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lossy materials, that accounts for both the nonlocal- 
bulk and local-surface nonlinear polarization sources. 
The problem is treated in the frequency domain under 
the undepleted-pump approximation. In particular, we 
demonstrate that the contribution of the nonlocal-bulk 
source to the SH electromagnetic field can be taken into 
account by introducing an equivalent surface magnetic 
current on the particle boundary, leading to a great sim- 
plification of the numerical problem. Unlike ref. [l3[ , the 
numerical formulation has been derived in the framework 
of the distribution theory in order to correctly account 
for the discontinuity of the discrete representation of the 
SH sources. 

The paper is organized as follows. Section 2 presents 
the mathematical statement of the problem along with 
the expression of both the nonlocal-bulk and the local- 
surface nonlinear polarization sources of the SH radia- 
tion. Then Section 3 is devoted to the SIE formulation 
and its numerical solution through the Galerkin method. 
Section 4 deals with the validation of the method, ac- 
complished through an extensive comparison with the 
full wave SHG Mie theory for nanoparticles with spheri- 
cal shapes proposed in ref. [Hj]. Finally, the SHG for a 
gold triangular prism is presented. 



II. STATEMENT OF THE PROBLEM 

In this section, we formulate the problem of SH genera- 
tion from a lossy NP of arbitrary shape in a homogeneous 
medium, when it is illuminated by a time-harmonic elec- 
tromagnetic field at frequency uj. The domain of the elec- 
tromagnetic field is the entire space R 3 , which is divided 

o 

into the interior of the metal domain Vi, the external 

o 

medium V e and the interface S. The closed surface S is 
oriented such that its normal n points outward. Further- 
more, r denotes the position vector with respect to an ar- 
bitrary reference point O. Since the intensities of the SH 
fields generated by noble metals NP are always order of 
magnitude weaker than the intensities of the pump fields, 
we can assume that the SH fields do not couple back 
to the fundamental fields {undepleted-pump approxima- 
tion). Therefore, the non-linear electromagnetic problem 
can be decomposed in two linear scattering problems at 
the fundamental frequency and at the SH frequency, and 
the nonlinear response of the material is phenomenolog- 
ically described by SH polarization sources induced by 
the fundamental electric field. 



We denote with E 



(e< w \h?">) 



the fundamental fields 



at frequency uj in the region Vi with I — i,e, and with 



f2 and the permeability of the metal and with and (e e , /i e ) 
the permeability and the permittivity of the embedding 
medium. In order to calculate the SH radiation gener- 
ated by the metal particle we have to evaluate: 1) the 
electric fields Ej at the fundamental frequency induced 
by the pump electromagnetic field; 2) the bulk and sur- 
face nonlinear sources generated by E^ and 3) the SH 
fields n3j 2w \ Hp 1 "^ generated by the nonlinear sources. 

Aiming at the solution of scattering problem at the 
fundamental frequency, we introduce the scattered fields 
w( w ' s ) ; fj| w,s ' ) with I = i,e, defined as: 



E 



E (w,s) _ E (w) 



, , in V e < . . in V> 

The scattered fields at the fundamental frequency can be 
determined by solving the problem: 



V x E^' s) = -jw W H, M 



in Vi, with I = i,e 



n x E 



(u>,«) _ F K«)A _ _ 



n x ( 



H 



E 

(ui,s) _ Jj( w > s ) 



= n x E 



n x H 



o 



onS 



(2) 

with the radiation condition at infinity. 

Once the solution of problem [2] has been obtained, the 
SH polarization sources can be calculated. Noble met- 
als are centrosymmetric materials, therefore the leading 
term of the bulk nonlinear polarization density P^ 2 "' is: 



,(2u0 



>(2) 



; 2) :E^fvE^ for ref, (3) 



where \b i s the quadrupolar contribution to the second- 
order nonlinear bulk susceptibility of the metal. Noble 
metals also are isotropic, in the case of homogeneous ma- 
terial the relation [3] becomes 



pf") = £o7 (2) V (e.^ ■ E< w) ) + s Q S^ (e 2 (w) • V) E^\ 

(4) 

where and 6^ are material parameters. Measure- 
ments indicate that the second term in Q] is negligible 
with respect to the second one [l6j]. Therefore, the bulk 
polarization source is: 



P<*0 = eo7 Wv (e< 



E 



(2oi) jj(2w) 



( 

convention a(r) = 3?{A^ n ) (r) exp (jflt)} for represent- 
ing a time harmonic electromagnetic field. Furthermore, 

we denote with ^E W \ H^^j the pump fields incident in 

o 

V e , with (si {fl} , fii) the linear permittivity at frequency 



the SH field at frequency 2u>. We use the The local-surface SH polarization source P^ is 



p(2u>) _ <+( 2 ) . F (w) F (a;) 

F s - £ oXs ■ E, E i , 



(5) 



(6) 



>(2) . 



where \s ^ s the second-order surface nonlinear suscep- 
tibility. Since the nanoparticle surface has an isotropic 
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symmetry with a mirror plane perpendicular to it, the 

<->(2) 

surface-susceptibility tensor \s nas onr y three non- 

(2) (2) 

vanishing and independent elements, Xj_j__l> X_M|i| anc ^ 

(2) 



(2) 

Am 



X || II j. i where _L and || refer to the orthogonal 
and tangential components to the nanoparticle surface. 
Therefore, we have: 



p(2w) _ , 



X 



nnn 



,(2) 



(ntiti + nt 2 t 2 ) + 



(2) 
ll-LII 



(tinti + t 2 nt 2 ) 



(w) T -,(o;) 



(7) 



where (ti, t 2 , n) is a system of three orthonormal vectors 
locally defined on the particle surface. The contribution 
of tangential and normal component of the surface non- 
linear polarization are taken into account by introduc- 
ing the surface electric and magnetic ftm^ current 
density : 



^=j2^\ 



J2u) 



SO 



V7 n(2") 



(8a) 
(8b) 



where j = n x n x P^ 2 "^ and P'g^ — P^ 2 "' • n. 

Once the phenomenological sources are known, we 
can determine the SH fields by solving non-homogeneous 
problem: 



V x E< 2w) 

V x Hi 2tj) 



V x H 

n x (e< 2 ") 
n x (Hi 2 ^ 



j2ue e E^ 
-j2uj^i 1 
(2w) -j2«e i {2w}Ef 



in V e , 



in Vi, 



■j2uP b 



E 



(2u) 



H 



(2w) 



(2") 



on5 



(9a) 
(9b) 

(9c) 



with the radiation condition at infinity. Due to the lin- 
earity of the problem I9"all9"cl its general solution is the 
sum of the solution of the homogeneous equation and one 
particular solution of the inhomogeneous equation. The 
equations [Ha] and [9b] are satisfied by a vanishing magnetic 
field and by the electric field: 



E. 



(2w) 
part 



o 



so7 <2) V ( 



E M . E (w) 



if r e V e , 
if reVi. 



(10) 



Therefore, the solution of the non-homogeneous problem 
defined by eqs. I9"all9"cl can be written as: 



E 



H 



(2") 
I 

(2") 



E 
H 



(2") 



E 



(2") 
part 



(2w) 



with I = 



i,e , 



(11) 



problem: 

J V x E ; (2w) = -j2w W H| 2w) 
1 VxH| 2w) = j2wej{2w}E, (2u) 



in Vi with Z = i,e 



(2u) _ -p(2M\ = __(2u>) _ _(2«) 
m 5 



n x f E^ J - E 
nx (H^-H^l =tt1 2 ^ 



with the radiation condition at infinity, where: 

(2ui) 1 _ n (2u) 



and 



(2 W ) _ £q7 (2) / w (u,) „(u>) 



( E M. E M) 



on S , 
(12) 

(13) 
(14) 



In conclusion, the nonlocal bulk source can be equiva- 
lently treated by using a superficial magnetic current, 
associated to a fictitious surface polarization PP^n. 



III. SURFACE INTEGRAL EQUATIONS 

Let us define, for the same geometry described in the 
previous section, an auxiliary scattering problem that 
includes both the fundamental and the SH problems 



as special cases. Let 7Tq ' , tv 



(n,e) (n,m) 



be the impressed 



electric and magnetic sources on S. We denote with 
I = i,e the fields solution of the prob- 



lcm: 



(f!,s) Ti(n,s) 



V x E 

V x H 



(fi,s) _ 

(n,«) 



(n, s ) 



n x E 



? <n,s) 



E 



nx H 



r (n, s ) 



H 



(f2,s) 



in V; with . 



j-^ q M 

+7T 



on S , 



. (15) 

with the radiation condition at infinity. We obtain the 
fundamental problem described of eq. [H by substituting 
Q = ui and 



^0 

(u;,m) 



n x H, 



(a,) 



J 



n x E, 



(a;) 

o • 



(16) 



Analogously, we can obtain the SH problem described by 
eq. [T2]by substituting fi = 2to and 



(2o,e) _ ( 2w ) 



(2a 

- v 



= ^(2*) 



r (2w) 



(17) 



where ( E 



(2u>) jj(2w) 



We now derive the surface integral formulation for the 
problem of eq. [T31 By invoking the Love's equivalent 

with I = i,e is the solution of the principle 17 1 for the exterior medium, the region Vi 
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is filled up with the same material of the region V e , 
and the sources are removed. The equivalent currents 

(je n,e \ je^'^j , positioned on the external page S e of the 
surface 5, produce the field ( E p\Hi as) ) in the region 

o o 

V e , and null fields in the region Vi, i.e.: 



if re 1/,, 



Ei n ' s) (r) ifrey e , 



^ ) {j^ e) j^ ) }w=r (n , s) 



if r € Sli, 
(r) ifref e , 



where the sources (j^'^,^'™^ are given by: 



( :(n,e) 

Je 

s(fi,m) 



n x H< n ' 8 > 
n x Ei°' s) 



(18) 



(19) 



An equivalent problem may be set up also for the region 

o o 

V i. In this case the region Vi is filled up with the ma- 
terial with electromagnetic parameters (e; {51} , fii). The 

equivalent currents (if 1 ^, j| n ' m ^ , defined on the inter- 
nal page Sj of the surface S, produce the original fields 

o o 

in the region Vi and null field in the region V e , i.e. 



f 5 (n,e) ;(n,m) 
°i \->i 'J* 



v (n) f.(n,e) .ca,m) 



}(r) 



Ef' s) (r) ifreV, 
if r e V e 

Hf 2 ' s) (r) ifref, 







if r € V P- 



where the sources 



(20) 



Ji J are given by: 



on the surface S we obtain: 



{jf' e) jf' m) }) 



S» 



_ n x j(n,m) 


= 0, 


(23a) 


+ n x $W 


= 0, 


(23b) 


+ nxjf' m > 


= 0, 


(23c) 


-nx 


= 0, 


(23d) 



where we have defined the operators: 



H 



' 4) {jW,jW} = -nxnxjf {j«J { "0} 



(24) 

Eqs. I23al I23cl are the T-electric field integral equations 
(T-EFIE), eqs. [23dl[23bl are the T-magnetic field integral 
equations (T-MFIE). By subtracting eq. I23cl from eq. I23al 
and eq. I23dl from I23bl and substituting the expressions 
of the operators , provided in the appendix [3] 

we obtain the PMCHWT formulation: 

C<n> x («) = y <n>, (25) 
where we have defined the operator 



CeTe^ 








_ K (n,t) 






>-i^-(n,i) 

M ' i 


1 





i 








1 





I 



c (n) = 



(26) 

the vector of unknowns x^- 1 , and the excitation vector 



Je 

• (Q,m) 
Je 
.(fi,e) 
J« 

.(n,m) 



1 (fi,e) 



7T. 



(fl,e) 


(n, ; 



, (27) 



( .(n,e) tt(£2.s) 

j =-nxH ; 



.(O.m) 



n x E 



(Q,s) 



(21) 



Furthermore, the sets of equivalent currents 

f.(Cl.e) .(n,m)\ i .{n.m)\ . -, 

(je ,Je J and Ijj- ,ji I are not indepen- 
dent. In fact, combining the boundary conditions of the 
problem [TBI and the definitions [T9l and |2"T1 we obtain: 



• .(O.e) .(n, e ) _ (n,e) 
Ji "T Je — 77 ' 

.(£2,m) , ; (n,m) _ _(n,m) 
. Ji "r Je — ^0 



(22) 



We now consider the limit of eqs. [TH1 and [201 a s the 
point r approaches the surface S from its external face S e 
and internal face Si. Projecting the resulting equations 



and Q {£1} = y/iu/ei {^}, with I = i,e. The expression 

of the operators 7^ (n '* ) and K-e can be found in the 
appendix lAl The expression of the excitation vector does 
not agree with the corresponding expression provided in 

m. 



IV. NUMERICAL SOLUTIONS 

In the present section, we describe the numerical meth- 
ods used to solve the PMCHWT formulation of eq. [25] 
both at the fundamental and at the SH frequency. In 
particular, we have exploited the Galerkin testing proce- 
dure to solve the singular SIEs, choosing the set of basis 
functions to coincide with the set of test functions [l7[ . 

Let us introduce a triangular mesh with N e edges. We 
denote with M. the mesh surface, the length of the p th 
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FIG. 1. (a) RWG basis function associated to the p t edge 
and denned on the triangle pair Tp ,Tp . In Tp and Tp f p 
is proportional to the vector (r — v^) and (v~ — r), respec- 
tively, (b) For the triangle Tp we show the support of the 
curve jp its tangent vector tp and binormal vector . 



edge as l p , and the two triangles sharing the p th edge as 
T+ and T~ (Fig. Q] (a)). They have an area of and 
A~ respectively. The vertices of T+ and T~ , which are 
opposite to the p th edge, are indicated by v+ and v~ re- 
spectively. The unknown electric and magnetic currents 
on M. are expanded in terms of the Rao- Wilton- Glisson 
(RWG) functions set f\, f 2 , . . . , fjv e . The RWG func- 
tion relative to the p th 



f P (r) 



fp-(r) 



edge is given by 




rer+ 

reT- 
otherwise 



(28) 



and the corresponding support is S p 
introduce also the symmetric product (f , g) = J s . 



T + U T~ 



We 
, f • gdS 

for the space of square integrable vector functions defined 
onS. 

We assembly the discrete operators associated to JCf 1 '^ 
and 7^ by projecting on the testing function f p the 

and 7^ n '*' of the ba- 



image through the operators ICf 1 '*' 
sis function f q . The above operators are explicitly derived 
in the section [B] of the appendix. Similarly, we assembly 
the discrete excitation vectors for both the fundamental 
and SH problem by projecting on the testing function f p 

the fiplrU 7r (w ' e) 7r ( "' m) 7r (2w) tt (2w) n x tt 

UK, I1C1US 7Tq , 7T n , 7T e , 7T m , 7T,, , II A 71, 



'0 



n X 7Tq , n X 7Te ' 



n x 7r 



(2") 



and n x 7r 



(2w) 



(2u>) 



( W ,e) 
' 

In par- 

(2w) 



ticular, in order to compute the quantities 7Te ', 7r^ 
7r[ 2w ' ) we need to evaluate the fundamental field on the 



internal page of the surface S : 



Ef > (r) 



e£? m - + JC ( r ) 



Si 



where 



E 



ft' W 
>) 



= -n x j 



.5, 



Si 



,- Vg-Ji ( r ) 

W£j {w} 

s («,e) 



(29) 



(30) 



piecewise constant function presenting a discontinuity of 
the first kind on the edges of the triangles. As a conse- 
quence, the calculation of the surface gradients V5P, 



(2w) 
S,n ' 



VsPg 2 ^ needed to obtain TTm"\ has to be carried 

out in the framework of the distribution theory. 

We now introduce the following notations. For each 
edge of the mesh, we define the two curves 7+ and 7" , 
having as support the perimeters of the triangles and 
, respectively, and are oriented counterclockwise when 
viewed from the end of the vector n; tp are their tangent 
vectors, and nip =tj xn are their binormal vectors (or- 
thogonal to both tp and n) (Fig. D](b)). The operators 
{•}| 7 and {'}| 7pi associate to a given physical quan- 
tity defined on A4, with the exception of the mesh edges, 
its limiting value as the evaluation point approaches the 
curve from the external or the internal side of the 
triangle T^, respectively. 

The gradient of Pg 2 ^ on the support of the basis func- 
tion fp, which is S p = T+ U T~ , has the following expres- 
sion: 



P 



(2u) 



P 



S.n 

(2u 



S,n 



7p + = 


p(2co) 
r S,n 




) <5 7+ m 




p (2u>) 
r S.n 




) <V m 


7p,e 


Tp ) i> 



if r e Tp U Tp , 



(31) 



where the functions 5 + and S - are Dirac function im- 

7p 7p 



pulsive on the support of the curves 7+ and 7" , respec- 
tively, and elsewhere on Ai, such that JJ M S^±dS = 1. 
It is worth noting that the first two branches of eq. [31] 
assume the same values along the edge l p shared by the 

two triangles T+ and T~ . The superficial current TTm^ 
on the support of f p will be: 



if r € 7+, 

if r e Tp U Tp , 

where the line currents ip^ are defined as: 



(32) 



</>+ = 


1 

£0 


f p(2w) 


7p + . 


p(2w) 
r S,n 








7p,e 


p(2u) 



t+ 



(33) 



In conclusion, as a consequence of the linearity of the 
RWGs, the surface current TTm^ degenerates into line 
Since the equivalent electric current (r) are repre- currents V™ defined on the mesh's edges, directed along 

sented through the linear RWGs functions, E^J (r) is a the edge itself. 
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Similarly, in order to compute the magnetic surface 
currents tt^^ introduced in eq. I13( we need to calcu- 
late the gradient of Pf; 2 ^ ■ P& 1S a piecewise quadratic 
function presenting a discontinuity of the first kind along 
the mesh edges. Therefore, we can repeat the derivation 
of the previous paragraph, obtaining: 



(2u>) _ 



if re 7+, 
if re 7-, 

o + o — 

if r e T p U T p 



(34) 



where 



4>t 



1 

£0 
1 



P 



(2u0 



P, 



S,fc 
(2") 



P. 



(2w) 



S,6 



7p,. 



7j>,. 



S,i> 

3 (2«) 
S,6 



(35) 



Thus, the current tt£ is a superposition of surface cur- 
rents defined on the interior of the mesh's triangles and of 
line currents ipf defined on the mesh's edges and directed 
along the edge itself. The line currents have not been 
considered in [l3| , as a results the provided expression of 

the discrete excitation vector associated to (fp^m ) is 
incorrect. In the appendix [Cl all the discrete excitation 
vectors are explicitly calculated. 



V. DISCUSSION 

In the first part of this section, we validate our for- 
mulation for the problem of the SH generation from an 
isolated gold sphere using the full wave SH Mie solution 
(SH-Mie) of the problems and QU The SH-Mie solu- 
tion has been obtained by expanding the pump field, the 
non-linear sources and the SH fields in scries of vector 
spherical wave functions (TBj . Then, we investigate the 
SH generation properties of a triangular nano-prism. The 
scatterers are homogeneous gold particles, embedded in 
air. The investigated physical quantity is the SH radiated 
power per unit solid angle collected along the direction k 



dP^/dVL (k) 



lim 

'||— 5- + 00 



2Ce m 



e£°> (k) 



(36) 



From both theoretical and experimental studies it rc- 

(2) 

suits that the component Xj_|||| om y weakly contributes 
to the SH generation from noble metals [TH, . Follow- 
ing Bachelier et al. [l9[ we express the parameters X±±±> 
XmjjI' an< ^ 7^ through the dimensionless phenomenolog- 
ical Rudnick- Stern parameters a, b and d [20j, as follows: 



(2) (2) (2) 





a 


b 


d 




4' 


2" 


8_ 



v (l) (, .\ u l e ° 

Xb l W i 2 ' 

oj z eno 



where x£ ( w ) = ( e * ( w ) / e o — 1) i s the bulk linear sus- 
ceptibility of the metal, e is the electron charge (abso- 
lute value), no is the number density of the conduction 
electrons, lu p is the plasma frequency and eo is the per- 
mittivity of the vacuum. By choosing a = 1, b = — 1 
and d = 1 we obtain the Rudnick-Stern hydrodynamic 
model [20| . The gold dispersion is modeled by fitting 
experimental data [2l| . The integrals involved in the nu- 
merical solution have been evaluated using four Gauss 
quadrature points for each triangle of M. , and their sin- 
gularities have been managed exploiting singularity ex- 
traction techniques (22j. 

The gold sphere is excited by a plane wave of unitary 
intensity, namely Eg = exp (— ju!^/e e fi e z^j x propagat- 
ing along the positive direction of the z-axis and linearly 
polarized along x. The exciting wavelength corresponds 
to the plasmonic resonance of a gold spherical particle in 
the Ray leigh limit, i.e. 520nm. The relative permittivity 
used in the calculations are &j (520nm) = —3.88 — 2.63j 
and Si (260nm) = —1.20 — 4.67j for the fundamental and 
the SH problem, respectively. The mesh has a number 
of edges equal to iV e = 3747. The maximum expansion 
order of vector spherical wave functions, used in the cal- 
culation of the SH-Mie solution, is N = 10. Aiming at an 
accurate validation of the proposed approach, we show in 
fig. [5] the quantity dP' 2 ") / dQ as a function of the incli- 
nation angle 9 for = and = 90. These quantities 
have been calculated by the SH-Mie (orange triangles for 
= and red circles for — 90) and by the SH-SIE 
(continuous black line for = and continuous blue line 
for = 90). Several sphere's diameters D have been 
considered: the first row of fig. [2] (panels a-d), is rela- 
tive to an electrically small particle with D = 20nm, the 
second row (e-h) to D = lOOnm, the third row (i-1) to 
D = 200nm. We considered in the first three columns 
only one SH source at a time, whereas in the fourth col- 
umn the different nonlinear polarization sources are si- 
multaneously active. In particular, in the first column 
of fig. [H corresponding to panels (a,e,i), we considered 
only the nonlocal-bulk source P fe „ using the coefficients 



(2) 

7>X mi 



p) \ _ 



( 

only the surface source P s t 



^ = (1, 0, 0). In the second column (b,f,j) 
' has been considered as- 
(0, 1, 0), whereas in the third 



summg ^7jX||j_||>X 

column (c,g,k) only the surface source P, 
considered assuming ^7 
tually, in panels (d,h,l,p 



(2u 



has been 



S,n 

.x[| 2 2n,xi 2 ix) = (0,0, i). Even " 

! we considered a realistic physi- 
cal case in which the different weights are calculated us- 
ing eq. [37] with (a,b,d) = (1, —1, 1). In all the different 
cases shown in fig. [5] very good agreement between the 
SH-SIE and the SH-Mie formulations can be found. 

In Fig. [3] we calculate the relative error £ of the SH- 
SIE formulation with respect to the SH-Mie formulation, 
defined as: 



(37) £ = (dP^/dn - p^/dn) I (p^/drij 



(38) 
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FIG. 2. SH radiated power per unit solid angle (expressed in Watt/Steradian) of a gold sphere with diameter (a-d) 20nm, (e-h) 
lOOnm (i-1) 200nm. Different combinations of SHG polarization sources have been considered, namely ^y' 2 ' , x[|x|| > Xxxx) = 
(1,0,0) in panels (a,e,i); (0,1,0) in panels (b,f,j); (0,0,1) in panels (c,g,k). In panels (d,h,l,p) the sources have been weighted 
using eq. (37jwith (a, b, d) = (1, —1, 1). The sphere is excited by a x-polarized plane wave of unitary intensity propagating along 
the positive z axis with wavelength 520nm. The diagrams are in linear scale. 




inclination angle 6 inclination angle 6 

FIG. 3. Relative error of dP {2uj) /dft calculated with the SH- 
SIE with respect to the SH-Mie solution for a sphere of di- 
ameter D = lOOnm as a function of the inclination angle 
9 for (a) 4> = and (b) <j> = 90. Different combinations 
of SHG polarization sources have been considered, namely 

(7 (2) ,X||5||,Xx 2 ix) = (1.0,0) (black line); (0,1,0) (red line); 
(0,0, 1) (blue line); the hydrodynamic model of eq. 1371 with 
(a, b, d) — (1, — 1, 1) (green line). 




FIG. 4. Radiation diagram of dP f2uj ^ jdQ, (expressed in 
Watt/Steradian) for a gold sphere with diameter: (a) 20nm, 
(b) lOOnm (c) 200nm. The sphere is excited by a x-polarized 
plane wave of unitary intensity propagating along the posi- 
tive direction of the z-axis with wavelength 520nm. The radi- 
ated powe is evaluated at 260nm. The diagrams are in linear 
scale. The sources have been weighted using eq. [37] with 
(a, b, d) = (1,-1,1). 



where Pfljg /dSl is the radiated power per unit solid an- 
gle calculate with the SH-Mie formulation. The error is 
plotted as a function of inclination angle 6 for (a) = 
and (b) cj> = 90 for of the for the D — lOOnm sphere. All 



the investigated cases present a comparable error, below 
3%. Only in correspondence of the forward and back- 
ward scattering direction, when the corresponding value 
of dP( 2uj ^ /d£l is very small the error is higher. 

In figure[3]we show the radiation diagram 
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FIG. 5. (a,c) Magnitude of the fundamental at the fun- 
damental frequency on the surface of the triangular nano- 
prism excited by a monochromatic (A = 690nm) plane wave 
of unitary intensity propagating along the positive z — axis 
and linearly polarized along (a) x and (c) y, and correspond- 
ing radiated power per unit solid angle (inset). (b,d) Radi- 
ation diagram of dP^ 2 "' /dQ, (expressed in Watt/Steradian) 
for both a (b) x-polarized (d) y-polarized pump. The sources 
have been weighted using eq. [STJwith (a, b, d) — (1, — 1, 1). 



for different diameters of the sphere, D — (a) 20nm (b) 
lOOnm (c) 200nm. The SH radiated power vanishes for 
the forward and backscattering directions for all the in- 
vestigated diameters, because the sphere is rotationally 
symmetric around the propagation direction. As we in- 
crease the sphere's diameter we notice the appearance of 
multiple oscillations corresponding to high order modes. 
Moreover, as the particle size increases, the inclination 
angle corresponding to the first maximum SH intensity 
approaches the forward scattering direction. This be- 
havior has been experimentally demonstrated for silver 
spherical particles in a recent work by Gonella et al. [23[ . 

We now consider a triangular nano-prism of height 
h = 40nm with rounded edges of length 200nm. The 
radius of the osculating cylinder (sphere) at each edge 
(corner) is 10nm. This shape is of extreme interest in 
nanoplasmonics as the building block of bow-tie nano- 
antennas, whose SH properties have been recently in- 
vestigated in [24j , We chose as the exciting wavelength 
A = 690nm, corresponding to the peak of the extinc- 
tion cross section over the visible range. We first an- 
alyze the scattering properties at the fundamental fre- 
quency. In figure [5] (a,c) we show the magnitude of the 
total electric field at the fundamental frequency on the 
surface of the triangular nano-prism when the incident 
plane wave is polarized along (a) the x axis and (c) the 
y axis. When the incident plane wave is x-polarized 



two hot spots arise at the two ends of the edge paral- 
lel to the incident polarization, whereas in correspon- 
dence of a y-polarized excitation the remaining corner 
experiences the highest field enhancement. The radiated 
power dP^' /dSl at the fundamental frequency can be 
attributed to an electric dipole oriented along the polar- 
ization direction (inset of Fig. [5](a,c)). In panels (b,d) 
we study the SH radiated power dP^- 2 ^ /dfl by the tri- 
angular nano-prism. The different sources are weighted 
using eq. [37] with (a, b, d) = (1, —1, 1). Unlike the linear 
scattering, dP' 2 "' / 'dVL shows a richer behavior dominated 
by higher order modes. Moreover, as the analyzed nano- 
prism is not rotationally symmetric around the propaga- 
tion direction, the radiated power in both the forward 
and backscattering directions is non zero. Furthermore, 
the shape of dP^ 2u) ^ / dH. is affected by the distribution 
of the electric field on the surface of the particle, shown 
in fig. [5] (a,b). In fact, we attribute the presence of the 
two lobes in panel to the strong field localization that oc- 
curs in correspondence of the two corners of the edge of 
the triangular-nanoprism parallel to x-axis. In contrast, 
when the incident field is directed along the y-axis, the ra- 
diation diagram features one lobe driven by the hot-spot 
on the remaining corner. Although both the surface area 
and the volume of the nano-prims are lower than those 
of the sphere of diameter 200nm investigated earlier, the 
maximum intensity of SH radiated power per unit solid 
angle dP^ 2 ^ / 'dQ of the nano-prism (Fig. EJb,d)) is al- 
most two order of magnitude higher compared to the 
sphere (Fig. S (b)). This fact is due to both the higher 
values of the non-linear susceptibilities of gold at 690nm 
and to the NP asymmetric shape 25 1. 



VI. CONCLUSIONS 

We have presented a surface integral equation method 
to accurately solve the second harmonic generation prob- 
lem in metal nanoparticles with arbitrary shapes ac- 
counting for both the nonlocal-bulk and local-surface 
nonlinear polarization sources. We have demonstrated 
that the contribution of the nonlocal-bulk nonlinear po- 
larization source to the SH electromagnetic field can be 
taken into account by introducing an equivalent surface 
magnetic current on the nanoparticle boundary. We com- 
pared the numerical solution of the SH surface integral 
equations to the analytical Mie solution for spheres of 
several radii, obtaining very good agreement for any in- 
vestigated case. The present treatment could also be 
easily extended to modeling sum-frequency and higher 
harmonic generation. The developed method paves the 
way to a better understanding of the process of SHG 
in arbitrarily shaped nanoparticles. This approach can 
have a high impact in the design of novel nanoplasmonic 
devices with enhanced SHG emission, including for in- 
stance sensor probing physical and chemical properties 
of material surfaces. 
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In the present section, we provide the expression of the 
electromagnetic field generated by a surface distribution 
of electric and magnetic currents (j^' 6 -* , j( n ' m )) defined 
on S and radiating in a homogeneous medium with con- 
stitutive parameters (ez{fi},^;) at frequency ft. The 
electromagnetic field radiated by both the electric and 
magnetic magnetic currents is: 

E< n > (r) = £l Q) {j (0 < e) ,j (0 < m) } (r) (Ala) 
H (fi) (r) = n (n) j.(n,e) )j( fi, ro ) j (r) (Alb) 

where: 



£ (n) |.(n, e)>j (n, ro ) } = (; {Q}T (n) |. ( o, e ) j + ^n) | j(n , m) j 



n {a) | j( n,e) J (n,m)j 



_£(!>) |.(0,e)j 



/C ; (0) {w} (r) = - / w (r') x Vg^ (r - r') dS' (A3a) 
Js 

TP {w} (r) = -jk, {fi} f g^ (r - r') w (r') dS'+ 



jh {fl} Js 



(o, 

V' 5 f 2) (r - r') V' s ■ w (r') dS" 
(A3b) 



+ < 



+ 



2jfe,{0} 



[V s -j (0 ' e) ] n- inxj( n > : 



if r i S, 
if r e Si, 



[V s .j( a ' e )] n+ |n x j( n ' m ) ifreS e 
if r ^ 5, 



■±nx j 



(o,e) ( r ) + &r,: ) , ( ^ n if row 



2j( l {n}k l {n} 

-in x i( f2 - e ) Crl - [ Vs -j ( "' m) W] n ;f rf S 
2 nxj yi) 2 jCi{n}k l {n} n 11 1 t o e 



(A2) 



I 

Vs- denotes the surface divergence, JCi and Ti denote obtain: 
the integral operators: 



Irt In 



4 r =± s =± ^P" 4 ? 



dS (Bl) 



is the homogeneous space Green's function, i.e. 



^(n.t) {.} = _ n x n x ^(n) |.| 
] {•} = -n x n x 7^ {•} 



(A4a) 
(A4b) 



Appendix B: Discrete Operators 



Similarly, by projecting on the testing function f p the 
image of the basis function f q through the operators 
J~( n -^ W e obtain: 



, f T (fi,t) n __, , \ \ ^ rs 

\Ip,J-l Iq) — tpt g / j 2_ J ArAs 
r=± s=± 



By projecting on the testing function f p the image +_ — jj jj ( r — r') dS'dS 



through the operators K\ ' of the basis function f g we 



(B2) 
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Appendix C: Excitation Vectors 



The discrete excitation vector associated to the super- 
ficial electric current density ■k^'^ is: 



(f p ,nX7if" 



fp • n x Tvf u) dS 



j2u 



r=± 



P^-nxf^S (CI) 



where we have split the integral domain S p into the two 
constitutive triangular facets I~. 

We now assembly the discrete excitation vector asso- 
ciated to the superficial magnetic current density 7Tm"': 



is. 



,, » /. ,., ■ (C2) 



(f p , n x 7r 

Splitting the integral domain S p into the two constitutive 
triangular facets T p and using the expression of tt 
given in eq. I32| we obtain 



(2") 



(fp,nx7r^)} = 

l - / f+ • n x ip+dl + i / f" • n x $- m dl (C3) 

Tp 7p 

where we have used the sampling property of S p , when 
they are applied to the discontinuous function f p (r), 
which is in T p , and elsewhere. Using the defini- 
tion of t/>^ given in eq. I33J we obtain: 



(fp,n x 7r 
1 



f+ ■ m 

9c- 7" I P P 
1 

2^ 



p(2w). 



p(2w), 
^S.n <y+. 



fp m p 



7p 



P, 



(2w) | 



P, 



(2w) | 



<2Z (C4) 



The quantities f^ • vanish on 7^ except on the com- 
mon edge l p , where they are equal to ±1. Therefore, we 
have: 



lp 

£0 



p(2w), _ p(2w). 



(C5) 



where we have exploited the fact that -Pg 2 ^L± = 

■Pff = ^fi^lx*' and P S,n\ T ± is the value of Unc- 
tion Pj 2 "' 1 in the interior of 

Then, we derive the expression for the discrete excita- 
tion vector associated to the superficial magnetic current 



density 7r 



(2a 



(fp,n X 7T 



(2w)\ 



fp (r) • n x nf u) dS (C6) 



The current Tr^ 2 "^ is a superposition of surface currents 
defined on the interior of the mesh's triangles and of line 
currents i/jf as exemplified in eq. |34] 

(fp, n x 7rf u) ) =I 1+ I 2 = ^J2 ( f f P r ' n x ^ldl+ 
4Eff fp r -V s P£ w) d5 (C7) 



£0 



The evaluation of the integral ii follows the same steps 
of (fp, n x 7im )• To evaluate the term / 2 , we first apply 



the vectorial identity f p • VsP fe „ ^ = V5 



P h „ ; Vs • fp", then the divergence theorem on the first of 
the two resulting terms, obtaining: 



p(2li)) n r 

r b,n l p 



1 

E0 



1 

£0 



E 

r=± 



P^Vs-^dS (C8) 



The first term on the r.h.s. of eq. IC8I is equal and oppo- 
site to I\ therefore they cancel in eq. IC71 i.e. 



fp, n x n^) = + -l p ^l- \ II P^dS (C9) 



£0 



r=± \ P 



,(2"), 



Now we now assembly the discrete excitation vector 



associated to the vector field n x tv 



(f p ,*™) = // t p -^dS = 



j2w E / / ? ' p s 2 t ^ (CIO) 



Then, we calculate the discrete excitation vector asso- 



ciated to the superficial vector field n x tv 



(2"). 



(fp,7r^))= // tp-n^dS (Cll) 

JJSp 

Splitting the integral domain S p into the two constitutive 
triangular facets , and using eq. [35] we have: 

(fp, =\i fp + ■ ^dl + fp" ■ ^ m dl (C12) 

** lp lp 

Using the eq. [33] we obtain: 



(fp^H^E/ p s:\j r p-%di 

Z£ o r=± J % 

-^E p £ ) k < / r ?-t;* (ci3) 



r=± 
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where we have used the fact that Pq 2 ^\ ± is constant 
along "fp. By noting that $ r • t£ dl = 0, we eventually 
obtain: 

(fp ' ^ = 2^ ^ f. Ps^kj; ■ % ^ (CM) 

^ r— ± 



The integral in eq. IC14I can be evaluated analytically. 

We now assembly the discrete excitation vector asso- 
ciated to n x 7r 



(2"). 



If _( 2liJ )\ 



tp ■ n^dS (C15) 



Splitting the integral domain S p into the two constitutive 
triangular facets and using eq. |33]we have: 



Z r=± J % 

4E// nxr p .VsP^dS (C16) 



£0 — 1 J JT r 

r=± j p 



By using eq. |35l the evaluation of the integral 1\ follows 
the same steps of (f p , 17m )'■ 



f r ■ t r 



P 



(2u>) 



b.n 



D (2u) 
6,n 



(C17) 

Concerning the term / 2 > we can apply the vectorial iden- 

tity [n x f£] • V S P^ ] = V s ■ [pg? (n x f;)] because 

the quantity V '$ • [n x f£] vanishes in 2J. Then, by us- 
ing the divergence theorem on the resulting integral, we 
obtain: 



h = 



1 E 



L p L p r b 



(2w) 



t/5 



(C18) 



Combining in eq. IC16I the derived expressions for the 
terms 1\ and I2 we obtain: 



(2a 



1 



2 EM? tt J)[Ck.+ i « 

r— ± "^p 



(2w), 
n Hi 



(C19) 
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